iSSF Vegetation Model Plotting

Here we are plotting the results of integrated step selection function (iSSF) models fitted to water buffalo in Arnhem Land. These models contain vegetation categories as a covariate, with interactions for movement parameters.

Load packages

library(tidyverse)
library(TwoStepCLogit) # to make population estimates
library(sjPlot) # to plot model estimates
library(glmmTMB)# for the mixed models 
library(tictoc)
library(beepr)
library(lemon)
library(circular)
library(amt)

Import data

ssf_dat <- read_csv("pheno_start_end_dat_ssf10rs_2024-03-14.csv")

New names:
Rows: 1062908 Columns: 42
── Column specification
──────────────────────────────────────────────────────── Delimiter: "," chr
(4): pheno_start, pheno_end, season, day_period dbl (35): ...1, id, burst_,
x1_, x2_, y1_, y2_, sl_, ta_, dt_, step_id_, sr... lgl (1): case_ dttm (2):
t1_, t2_
ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
Specify the column types or set `show_col_types = FALSE` to quiet this message.
• `` -> `...1`

# ssf_dat <- read_csv("pheno_forest_start_end_dat_ssf10rs_2025-03-02.csv")

Prepare the data

# scale the covariates
ssf_dat <- ssf_dat %>% mutate(srtm_start_scaled = scale(srtm_start),
                            srtm_end_scaled = scale(srtm_end),
                            water_dist_start_scaled = scale(water_dist_start),
                            water_dist_end_scaled = scale(water_dist_end),
                            ndvi_start_scaled = scale(ndvi_start),
                            ndvi_end_scaled = scale(ndvi_end),
                            nbr_start_scaled = scale(nbr_start),
                            nbr_end_scaled = scale(nbr_end),
                            sl_scaled = scale(sl_),
                            log_sl_scaled = scale(log_sl_),
                            cos_ta_scaled = scale(cos_ta_)
)

# get scaling factors
sl_scale_sd <- attr(ssf_dat$sl_scaled, "scaled:scale")
log_sl_scale_sd <- attr(ssf_dat$log_sl_scaled, "scaled:scale")
cos_ta_scale_sd <- attr(ssf_dat$cos_ta_scaled, "scaled:scale")

# rename the pheno categories to set a different reference level
ssf_dat$pheno_end <- relevel(as.factor(ssf_dat$pheno_end), ref = "shrub_grass")
ssf_dat$pheno_start <- relevel(as.factor(ssf_dat$pheno_start), ref = "shrub_grass")

unique(ssf_dat$pheno_start)

[1] shrub_grass   dry_grassland tree_grass    bare_ground   open_forest  
Levels: shrub_grass bare_ground dry_grassland open_forest tree_grass

unique(ssf_dat$pheno_end)

[1] shrub_grass   dry_grassland tree_grass    bare_ground   open_forest  
Levels: shrub_grass bare_ground dry_grassland open_forest tree_grass

ssf_dat <- ssf_dat %>% filter(!pheno_start == "saltpan_mudlfat")

# separate the buffalo by season 
ssfdat_wet <- ssf_dat %>% filter(season=="wet") # wet season
ssfdat_dry <- ssf_dat %>% filter(season=="dry") # dry season only 

Plot the used-available distributions

Dry season

# checking for all individuals
buffalo_ids <- unique(ssf_dat$id)

# inspect the data for dry season 
ssfdat_dry %>% 
  group_by(case_, pheno_end) %>% 
  summarize(n = n()) %>% 
  mutate(prop = n / sum(n), 
         label = paste0(round(prop * 100, 1), "%")) %>% 
  ggplot(aes(pheno_end, prop, fill = case_, group=case_,label = label)) + 
  geom_col(position = position_dodge2()) +
  geom_text(size = 4, vjust = -0.25, position = position_dodge(width = 1)) +
  labs(x = "Land use class", y = "Proportion", fill = "Case")+
  ggtitle("Dry season", subtitle = "All buffalo")+
  scale_fill_brewer(palette = "Paired", name="Case", 
                    breaks=c("FALSE", "TRUE"), labels=c("Available", "Used")) +
  theme_light() +
  theme(axis.text.x = element_text(angle = 20, hjust = 1))

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

ggsave("outputs/plots/used-available_all_dry.png", device = "png", 
       width = 150, height = 90, units = "mm", scale = 1.5, dpi = 1000)

Inspect the data for each buffalo

for(i in 1:length(buffalo_ids)) {
  print(ssfdat_dry %>% filter(id == buffalo_ids[i]) %>%
          group_by(case_, pheno_end) %>% 
          summarize(n = n()) %>% 
          mutate(prop = n / sum(n), 
                 label = paste0(round(prop * 100, 1), "%")) %>% 
          ggplot(aes(pheno_end, prop, fill = case_, group=case_,label = label)) + 
          geom_col(position = position_dodge2()) +
          geom_text(size = 4, vjust = -0.25, position = position_dodge(width = 1)) +
          labs(x = "Land use class", y = "Proportion", fill = "case_")+
          ggtitle("Dry season", subtitle = paste0("Buffalo ", buffalo_ids[i])) +
          scale_fill_brewer(palette = "Paired", name="Case", 
                            breaks=c("FALSE", "TRUE"), labels=c("Available", "Used")) +
          theme_light() +
          theme(axis.text.x = element_text(angle = 20, hjust = 1))
  )
}

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

Wet season

ssfdat_wet %>% 
  group_by(case_, pheno_end) %>% 
  summarize(n = n()) %>% 
  mutate(prop = n / sum(n), 
         label = paste0(round(prop * 100, 1), "%")) %>% 
  ggplot(aes(pheno_end, prop, fill = case_, group=case_,label = label)) + 
  geom_col(position = position_dodge2()) +
  geom_text(size = 4, vjust = -0.25, position = position_dodge(width = 1)) +
  labs(x = "Land use class", y = "Proportion", fill = "case_")+
  ggtitle("Wet season", subtitle = "All buffalo")+
  scale_fill_brewer(palette = "Paired", name="Case", 
                    breaks=c("FALSE", "TRUE"), labels=c("Available", "Used")) +
  theme_light() +
  theme(axis.text.x = element_text(angle = 20, hjust = 1))

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

ggsave("outputs/plots/used-available_all_wet.png", device = "png", 
       width = 150, height = 90, units = "mm", scale = 1.5, dpi = 1000)

For each individual

for(i in 1:length(buffalo_ids)) {
  print(ssfdat_wet %>% filter(id == buffalo_ids[i]) %>%
          group_by(case_, pheno_end) %>% 
          summarize(n = n()) %>% 
          mutate(prop = n / sum(n), 
                 label = paste0(round(prop * 100, 1), "%")) %>% 
          ggplot(aes(pheno_end, prop, fill = case_, group=case_,label = label)) + 
          geom_col(position = position_dodge2()) +
          geom_text(size = 4, vjust = -0.25, position = position_dodge(width = 1)) +
          labs(x = "Land use class", y = "Proportion", fill = "case_")+
          ggtitle("Wet season", subtitle = paste0("Buffalo ", buffalo_ids[i]))+
          scale_fill_brewer(palette = "Paired", name="Case", 
                            breaks=c("FALSE", "TRUE"), labels=c("Available", "Used")) +
          theme_light() +
          theme(axis.text.x = element_text(angle = 20, hjust = 1))
  )
}

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

`summarise()` has grouped output by 'case_'. You can override using the
`.groups` argument.

Read in the saved model objects

buffalo.pheno.dry.mvmt <- readRDS("buffalo_pheno_pooled_mvmt_10rs_shrubgrassref_dry_2025-05-20.rds")
buffalo.pheno.wet.mvmt <- readRDS("buffalo_pheno_pooled_mvmt_10rs_shrubgrassref_wet_2025-05-20.rds")

Model summaries

# Summaries are also printed below
# summary(buffalo.pheno.dry.mvmt)
# summary(buffalo.pheno.wet.mvmt)

Model fitting - interacting with movement

Dry season

tic()
buffalo.tmp <- glmmTMB(case_ ~ 
                         pheno_end + 
                         sl_scaled +
                         log_sl_scaled +
                         cos_ta_scaled +
                         pheno_start:(sl_scaled + log_sl_scaled + cos_ta_scaled) +
                         (1 | step_id),
                       family = poisson(), 
                       data = ssfdat_dry, 
                       doFit=FALSE,
                       control = glmmTMBControl(optimizer = optim, optArgs = list(method="BFGS")))

#' Set variance of random intercept to 10^6
buffalo.tmp$parameters$theta[1] <- log(1e3)
nvarparm<-length(buffalo.tmp$parameters$theta)
buffalo.tmp$mapArg <- list(theta=factor(c(NA)))

# to render the script we comment out the model fitting, as we are loading a saved model
# buffalo.pheno.dry.mvmt <- glmmTMB:::fitTMB(buffalo.tmp)

summary(buffalo.pheno.dry.mvmt) 

 Family: poisson  ( log )
Formula:          
case_ ~ pheno_end + sl_scaled + log_sl_scaled + cos_ta_scaled +  
    pheno_start:(sl_scaled + log_sl_scaled + cos_ta_scaled) +  
    (1 | step_id)
Data: ssfdat_dry

      AIC       BIC    logLik -2*log(L)  df.resid 
  1161944   1162171   -580952   1161904    620303 

Random effects:

Conditional model:
 Groups  Name        Variance Std.Dev.
 step_id (Intercept) 1e+06    1000    
Number of obs: 620323, groups:  step_id, 56393

Conditional model:
                                        Estimate Std. Error z value Pr(>|z|)
(Intercept)                            -2.943079   4.211029  -0.699 0.484616
pheno_endbare_ground                    0.207230   0.016675  12.428  < 2e-16
pheno_enddry_grassland                 -0.004613   0.016601  -0.278 0.781098
pheno_endopen_forest                   -0.028587   0.020086  -1.423 0.154663
pheno_endtree_grass                    -0.082072   0.018900  -4.342 1.41e-05
sl_scaled                               0.017588   0.008004   2.197 0.027998
log_sl_scaled                          -0.009450   0.008218  -1.150 0.250170
cos_ta_scaled                          -0.001376   0.006581  -0.209 0.834424
sl_scaled:pheno_startbare_ground       -0.032633   0.017720  -1.842 0.065540
sl_scaled:pheno_startdry_grassland      0.054589   0.016693   3.270 0.001075
sl_scaled:pheno_startopen_forest       -0.037419   0.017003  -2.201 0.027760
sl_scaled:pheno_starttree_grass         0.010700   0.018608   0.575 0.565258
log_sl_scaled:pheno_startbare_ground   -0.041416   0.015747  -2.630 0.008537
log_sl_scaled:pheno_startdry_grassland -0.095950   0.017085  -5.616 1.95e-08
log_sl_scaled:pheno_startopen_forest    0.058683   0.017622   3.330 0.000868
log_sl_scaled:pheno_starttree_grass     0.006182   0.019581   0.316 0.752240
cos_ta_scaled:pheno_startbare_ground   -0.046250   0.012817  -3.608 0.000308
cos_ta_scaled:pheno_startdry_grassland  0.009413   0.013829   0.681 0.496068
cos_ta_scaled:pheno_startopen_forest   -0.036834   0.013465  -2.735 0.006229
cos_ta_scaled:pheno_starttree_grass    -0.037956   0.014868  -2.553 0.010685
                                          
(Intercept)                               
pheno_endbare_ground                   ***
pheno_enddry_grassland                    
pheno_endopen_forest                      
pheno_endtree_grass                    ***
sl_scaled                              *  
log_sl_scaled                             
cos_ta_scaled                             
sl_scaled:pheno_startbare_ground       .  
sl_scaled:pheno_startdry_grassland     ** 
sl_scaled:pheno_startopen_forest       *  
sl_scaled:pheno_starttree_grass           
log_sl_scaled:pheno_startbare_ground   ** 
log_sl_scaled:pheno_startdry_grassland ***
log_sl_scaled:pheno_startopen_forest   ***
log_sl_scaled:pheno_starttree_grass       
cos_ta_scaled:pheno_startbare_ground   ***
cos_ta_scaled:pheno_startdry_grassland    
cos_ta_scaled:pheno_startopen_forest   ** 
cos_ta_scaled:pheno_starttree_grass    *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

toc()

3.878 sec elapsed

beep(sound = 2)

# plot result as incidence rate ratios
# plot_model(buffalo.pheno.dry.mvmt)
plot_model(buffalo.pheno.dry.mvmt, transform = NULL)

# saveRDS(buffalo.pheno.dry.mvmt, file = paste0("buffalo_pheno_pooled_mvmt_10rs_shrubgrassref_dry_", Sys.Date(), ".rds"))
# paste0("Model location: buffalo_pheno_pooled_mvmt_10rs_shrubgrassref_dry_", Sys.Date(), ".rds")

Wet season

tic()
buffalo.tmp <- glmmTMB(case_ ~ 
                         pheno_end + 
                         sl_scaled +
                         log_sl_scaled +
                         cos_ta_scaled +
                         pheno_start:(sl_scaled + log_sl_scaled + cos_ta_scaled) +
                         (1 | step_id),
                       family=poisson(), 
                       data = ssfdat_wet, 
                       doFit=FALSE,
                       control = glmmTMBControl(optimizer = optim, optArgs = list(method="BFGS")))

#' Set variance of random intercept to 10^6
buffalo.tmp$parameters$theta[1] <- log(1e3)
nvarparm<-length(buffalo.tmp$parameters$theta)
buffalo.tmp$mapArg <- list(theta=factor(c(NA)))

# to render the script we comment out the model fitting, as we are loading a saved model
# buffalo.pheno.wet.mvmt <- glmmTMB:::fitTMB(buffalo.tmp)

summary(buffalo.pheno.wet.mvmt)

 Family: poisson  ( log )
Formula:          
case_ ~ pheno_end + sl_scaled + log_sl_scaled + cos_ta_scaled +  
    pheno_start:(sl_scaled + log_sl_scaled + cos_ta_scaled) +  
    (1 | step_id)
Data: ssfdat_wet

      AIC       BIC    logLik -2*log(L)  df.resid 
 827723.7  827943.8 -413841.9  827683.7    442565 

Random effects:

Conditional model:
 Groups  Name        Variance Std.Dev.
 step_id (Intercept) 1e+06    1000    
Number of obs: 442585, groups:  step_id, 40235

Conditional model:
                                        Estimate Std. Error z value Pr(>|z|)
(Intercept)                            -3.078387   4.985386  -0.617   0.5369
pheno_endbare_ground                    0.524979   0.018049  29.086  < 2e-16
pheno_enddry_grassland                  0.122161   0.020143   6.065 1.32e-09
pheno_endopen_forest                    0.188558   0.025055   7.526 5.24e-14
pheno_endtree_grass                     0.056589   0.022771   2.485   0.0130
sl_scaled                              -0.002532   0.008825  -0.287   0.7742
log_sl_scaled                           0.150044   0.011151  13.455  < 2e-16
cos_ta_scaled                           0.084601   0.007900  10.709  < 2e-16
sl_scaled:pheno_startbare_ground       -0.021958   0.018241  -1.204   0.2287
sl_scaled:pheno_startdry_grassland      0.006751   0.021412   0.315   0.7525
sl_scaled:pheno_startopen_forest       -0.091184   0.022869  -3.987 6.68e-05
sl_scaled:pheno_starttree_grass        -0.040299   0.025016  -1.611   0.1072
log_sl_scaled:pheno_startbare_ground   -0.085673   0.018978  -4.514 6.35e-06
log_sl_scaled:pheno_startdry_grassland -0.050208   0.025023  -2.006   0.0448
log_sl_scaled:pheno_startopen_forest   -0.123207   0.022896  -5.381 7.40e-08
log_sl_scaled:pheno_starttree_grass    -0.160890   0.024680  -6.519 7.07e-11
cos_ta_scaled:pheno_startbare_ground   -0.105115   0.013989  -7.514 5.73e-14
cos_ta_scaled:pheno_startdry_grassland -0.014564   0.017465  -0.834   0.4043
cos_ta_scaled:pheno_startopen_forest   -0.159377   0.017228  -9.251  < 2e-16
cos_ta_scaled:pheno_starttree_grass    -0.142022   0.018050  -7.868 3.60e-15
                                          
(Intercept)                               
pheno_endbare_ground                   ***
pheno_enddry_grassland                 ***
pheno_endopen_forest                   ***
pheno_endtree_grass                    *  
sl_scaled                                 
log_sl_scaled                          ***
cos_ta_scaled                          ***
sl_scaled:pheno_startbare_ground          
sl_scaled:pheno_startdry_grassland        
sl_scaled:pheno_startopen_forest       ***
sl_scaled:pheno_starttree_grass           
log_sl_scaled:pheno_startbare_ground   ***
log_sl_scaled:pheno_startdry_grassland *  
log_sl_scaled:pheno_startopen_forest   ***
log_sl_scaled:pheno_starttree_grass    ***
cos_ta_scaled:pheno_startbare_ground   ***
cos_ta_scaled:pheno_startdry_grassland    
cos_ta_scaled:pheno_startopen_forest   ***
cos_ta_scaled:pheno_starttree_grass    ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

toc()

2.493 sec elapsed

beep(sound = 2)

# plot result as incidence rate ratios
# plot_model(buffalo.pheno.wet.mvmt)
plot_model(buffalo.pheno.wet.mvmt, transform = NULL)

# saveRDS(buffalo.pheno.wet.mvmt, file = paste0("buffalo_pheno_pooled_mvmt_10rs_shrubgrassref_wet_", Sys.Date(), ".rds"))
# paste0("Model location: buffalo_pheno_pooled_mvmt_10rs_shrubgrassref_wet_", Sys.Date(), ".rds")

Plot both models

# plot_model(buffalo.pheno.dry.mvmt)
# plot_model(buffalo.pheno.wet.mvmt)
# 
# plot_model(buffalo.pheno.dry.mvmt, transform = NULL)
# plot_model(buffalo.pheno.wet.mvmt, transform = NULL)

Create manuscript figures

Dry season

coef_df <- data.frame("model_covariate" = names(fixef(buffalo.pheno.dry.mvmt)$cond),
                      "Covariate" = c("Intercept", "Floodplain", "DryGrassland", "OpenForest", "TreeGrass", 
                                      "sl", "log(sl)", "cos(ta)", 
                                      "sl:Floodplain", "sl:DryGrassland", "sl:OpenForest", "sl:TreeGrass",
                                      "log(sl):Floodplain", "log(sl):DryGrassland", "log(sl):OpenForest", "log(sl):TreeGrass",
                                      "cos(ta):Floodplain", "cos(ta):DryGrassland", "cos(ta):OpenForest", "cos(ta):TreeGrass"
                                      ),
                      "Estimate" = coef(summary(buffalo.pheno.dry.mvmt))$cond[, "Estimate"],
                      "SE" = coef(summary(buffalo.pheno.dry.mvmt))$cond[, "Std. Error"]
) %>%
  
  mutate(
    # Calculate confidence intervals
    LCI_90 = Estimate - qnorm(1 - (1 - 0.90) / 2) * SE, 
    UCI_90 = Estimate + qnorm(1 - (1 - 0.90) / 2) * SE,
    LCI_95 = Estimate - qnorm(1 - (1 - 0.95) / 2) * SE,
    UCI_95 = Estimate + qnorm(1 - (1 - 0.95) / 2) * SE,
    LCI_99 = Estimate - qnorm(1 - (1 - 0.99) / 2) * SE,
    UCI_99 = Estimate + qnorm(1 - (1 - 0.99) / 2) * SE,
    # Compute p-values
    pvalue = 2 * pnorm(-abs(Estimate) / SE)
  )

# Add stars indicating the significance
coef_df$Significance <- sapply(1:nrow(coef_df), function(x){
  if (coef_df$pvalue[x] <= 0.001){
    return("***")
  } else if (coef_df$pvalue[x] <= 0.01) {
    return("**")
  } else if (coef_df$pvalue[x] <= 0.05) {
    return("*")
  }
})

# Remove the intercept term
coef_df <- coef_df %>% filter(Covariate != "Intercept")

# Add a column indicating the preference
coef_df$Preference <- ifelse(coef_df$Estimate > 0, "Preferred", "Avoided")
coef_df$Preference <- factor(coef_df$Preference, levels = c("Avoided", "Preferred"))

# Specify the order in which the coefficients should be plotted
order <- c(
  "Floodplain",
  "DryGrassland",
  "OpenForest",
  "TreeGrass",
  "sl",
  "log(sl)",
  "cos(ta)",
  "sl:Floodplain",
  "sl:DryGrassland",
  "sl:OpenForest",
  "sl:TreeGrass",
  "log(sl):Floodplain",
  "log(sl):DryGrassland",
  "log(sl):OpenForest",
  "log(sl):TreeGrass",
  "cos(ta):Floodplain",
  "cos(ta):DryGrassland",
  "cos(ta):OpenForest",
  "cos(ta):TreeGrass"
)

# Specify the order in which the coefficients should be plotted
order_no_mvmt <- c(
  "Floodplain",
  "DryGrassland",
  "OpenForest",
  "TreeGrass"
)

# Prepare dataset for plotting confidence intervals
coef_df2 <- coef_df %>%
  dplyr::select(Covariate, Estimate, Preference, LCI_90:UCI_99) %>%
  gather(key = confidence_level, value = value, LCI_90:UCI_99) %>%
  separate(col = confidence_level, into = c("Type", "Level"), sep = "_") %>%
  spread(key = Type, value = value) %>%
  mutate(Level = paste0(Level, "%"))

coef_df2_no_mvmt <- coef_df %>% filter(str_detect(model_covariate, "pheno_end")) %>%
  dplyr::select(Covariate, Estimate, Preference, LCI_90:UCI_99) %>%
  gather(key = confidence_level, value = value, LCI_90:UCI_99) %>%
  separate(col = confidence_level, into = c("Type", "Level"), sep = "_") %>%
  spread(key = Type, value = value) %>%
  mutate(Level = paste0(Level, "%"))

Plot the coefficients

min_x <- -0.4
max_x <- 0.65

# coefficients on the x-axis
dry_mvmt_plot <- ggplot(data = coef_df, 
                        aes(y = Covariate, x = Estimate, col = factor(Preference))) +
  geom_vline(xintercept = 0, color = "black", lty = 2, lwd = 0.3) +
  annotate(geom = "segment",
           x      = min_x, xend = max_x,
           y      = 12.5, yend   = 12.5,
           colour = "gray80", lty = 1, lwd = 0.3) +
  annotate(geom = "segment",
           x      = min_x, xend = max_x,
           y      = 15.5, yend   = 15.5,
           colour = "gray80", lty = 1, lwd = 0.3) +
  geom_point(shape = 3, size = 2.5) +
  # drop the filter argument when plotting all covariates
  geom_errorbarh(data = coef_df2,
                 aes(xmin = LCI, xmax = UCI, linewidth = factor(Level)), 
                 height = 0, alpha  = 0.5) +
  geom_text(aes(label = Significance, hjust = 0.5, vjust = 0), 
            show.legend = F) +
  scale_y_discrete(limits = rev(order)) + 
  theme_classic() +
  scale_x_continuous(limits = c(min_x, max_x), breaks = seq(min_x, max_x, 0.1)) +
  coord_capped_cart(left = "both", bottom = "both") +
  labs(x = expression(beta*"-Estimate")) +
  scale_color_manual(name   = "Preference", 
                     values = c("#FF8123", "#9B4200")) +
  scale_linewidth_manual(name   = "Confidence Level", 
                         values = c(2, 1, 0.3)) +
  # ggtitle("Dry season") +
  theme(legend.position   = "bottom",
        legend.margin     = margin(0, 50, 0, -20),
        legend.box.margin = margin(-5, -10, -5, -10),
        legend.text       = element_text(face = 3),
        legend.title      = element_text(face = 3)) +
  guides(colour = guide_legend(title.position = "top", title.hjust = 0.5),
         linewidth   = guide_legend(title.position = "top", title.hjust = 0.5, 
                                    override.aes = list(colour = "#FF8123")))

dry_mvmt_plot

Warning: Removed 8 rows containing missing values or values outside the scale range
(`geom_text()`).

ggsave(paste0("outputs/plots/pheno_model_dry_pooled_mvmt_", Sys.Date(), ".png"), device = "png", 
       width = 150, height = 150, units = "mm", scale = 1, dpi = 600)

Warning: Removed 8 rows containing missing values or values outside the scale range
(`geom_text()`).

Drop the filter argument when plotting all covariates

dry_mvmt_plot <- ggplot(data = coef_df %>% filter(str_detect(model_covariate, "pheno_end")), 
                        aes(y = Covariate, x = Estimate, col = factor(Preference))) +
  geom_vline(xintercept = 0, color = "black", lty = 2, lwd = 0.3) +
  geom_point(shape = 3, size = 2.5) +
  # drop the filter argument when plotting all covariates
  geom_errorbarh(data = coef_df2_no_mvmt,
                 aes(xmin = LCI, xmax = UCI, linewidth = factor(Level)), 
                 height = 0, alpha  = 0.5) +
  geom_text(aes(label = Significance, hjust = 0.5, vjust = 0), 
            show.legend = F) +
  # scale_y_discrete(limits = rev(order)) + 
  scale_y_discrete(limits = rev(order_no_mvmt)) + 
  theme_classic() +
  scale_x_continuous(limits = c(min_x, max_x), breaks = seq(min_x, max_x, 0.1)) +
  coord_capped_cart(left = "both", bottom = "both") +
  labs(x = expression(beta*"-Estimate")) +
  scale_color_manual(name   = "Preference", 
                     values = c("#FF8123", "#9B4200")) +
  scale_linewidth_manual(name   = "Confidence Level", 
                         values = c(2, 1, 0.3)) +
  theme(legend.position   = "bottom",
        legend.margin     = margin(0, 50, 0, -20),
        legend.box.margin = margin(-5, -10, -5, -10),
        legend.text       = element_text(face = 3),
        legend.title      = element_text(face = 3)) +
  guides(colour = guide_legend(title.position = "top", title.hjust = 0.5),
         linewidth   = guide_legend(title.position = "top", title.hjust = 0.5, 
                                    override.aes = list(colour = "#FF8123")))

dry_mvmt_plot

Warning: Removed 2 rows containing missing values or values outside the scale range
(`geom_text()`).

ggsave(paste0("outputs/plots/pheno_model_dry_pooled_mvmt_WITHOUT_MVMT_", Sys.Date(), ".png"), device = "png", 
       width = 150, height = 60, units = "mm", scale = 1, dpi = 600)

Warning: Removed 2 rows containing missing values or values outside the scale range
(`geom_text()`).

Wet season

coef_df <- data.frame("model_covariate" = names(fixef(buffalo.pheno.wet.mvmt)$cond),
                      "Covariate" = c("Intercept", "Floodplain", "DryGrassland", "OpenForest", "TreeGrass", 
                                      "sl", "log(sl)", "cos(ta)", 
                                      "sl:Floodplain", "sl:DryGrassland", "sl:OpenForest", "sl:TreeGrass",
                                      "log(sl):Floodplain", "log(sl):DryGrassland", "log(sl):OpenForest", "log(sl):TreeGrass",
                                      "cos(ta):Floodplain", "cos(ta):DryGrassland", "cos(ta):OpenForest", "cos(ta):TreeGrass"
                                      ),
                      "Estimate" = coef(summary(buffalo.pheno.wet.mvmt))$cond[, "Estimate"],
                      "SE" = coef(summary(buffalo.pheno.wet.mvmt))$cond[, "Std. Error"]
) %>%
  
  mutate(
    # Calculate confidence intervals
    LCI_90 = Estimate - qnorm(1 - (1 - 0.90) / 2) * SE, 
    UCI_90 = Estimate + qnorm(1 - (1 - 0.90) / 2) * SE,
    LCI_95 = Estimate - qnorm(1 - (1 - 0.95) / 2) * SE,
    UCI_95 = Estimate + qnorm(1 - (1 - 0.95) / 2) * SE,
    LCI_99 = Estimate - qnorm(1 - (1 - 0.99) / 2) * SE,
    UCI_99 = Estimate + qnorm(1 - (1 - 0.99) / 2) * SE,
    # Compute p-values
    pvalue = 2 * pnorm(-abs(Estimate) / SE)
  )

# Add stars indicating the significance
coef_df$Significance <- sapply(1:nrow(coef_df), function(x){
  if (coef_df$pvalue[x] <= 0.001){
    return("***")
  } else if (coef_df$pvalue[x] <= 0.01) {
    return("**")
  } else if (coef_df$pvalue[x] <= 0.05) {
    return("*")
  }
})

# Remove the intercept term
coef_df <- coef_df %>% filter(Covariate != "Intercept")

# Add a column indicating the preference
coef_df$Preference <- ifelse(coef_df$Estimate > 0, "Preferred", "Avoided")
coef_df$Preference <- factor(coef_df$Preference, levels = c("Avoided", "Preferred"))

# Prepare dataset for plotting confidence intervals
coef_df2 <- coef_df %>%
  dplyr::select(Covariate, Estimate, Preference, LCI_90:UCI_99) %>%
  gather(key = confidence_level, value = value, LCI_90:UCI_99) %>%
  separate(col = confidence_level, into = c("Type", "Level"), sep = "_") %>%
  spread(key = Type, value = value) %>%
  mutate(Level = paste0(Level, "%"))

coef_df2_no_mvmt <- coef_df %>% filter(str_detect(model_covariate, "pheno_end")) %>%
  dplyr::select(Covariate, Estimate, Preference, LCI_90:UCI_99) %>%
  gather(key = confidence_level, value = value, LCI_90:UCI_99) %>%
  separate(col = confidence_level, into = c("Type", "Level"), sep = "_") %>%
  spread(key = Type, value = value) %>%
  mutate(Level = paste0(Level, "%"))

Plot the coefficients

# coefficients on the x-axis
wet_mvmt_plot <- ggplot(data = coef_df, 
                        aes(y = Covariate, x = Estimate, col = factor(Preference))) +
  geom_vline(xintercept = 0, color = "black", lty = 2, lwd = 0.3) +
  annotate(geom = "segment",
           x      = min_x, xend = max_x,
           y      = 12.5, yend   = 12.5,
           colour = "gray80", lty = 1, lwd = 0.3) +
  annotate(geom = "segment",
           x      = min_x, xend = max_x,
           y      = 15.5, yend   = 15.5,
           colour = "gray80", lty = 1, lwd = 0.3) +
  geom_point(shape = 3, size = 2.5) +
  geom_point(shape = 3, size = 2.5) +
  geom_errorbarh(data = coef_df2,
                 aes(xmin = LCI, xmax = UCI, linewidth = factor(Level)), 
                 height = 0, alpha  = 0.5) +
  geom_text(aes(label = Significance, hjust = 0.5, vjust = 0), 
            show.legend = F) +
  scale_y_discrete(limits = rev(order)) +
  theme_classic() +
  scale_x_continuous(limits = c(min_x, max_x), breaks = seq(min_x, max_x, 0.1)) +
  coord_capped_cart(left = "both", bottom = "both") +
  labs(x = expression(beta*"-Estimate")) +
  scale_color_manual(name   = "Preference", 
                     values = c("#238DFF", "#004691")) +
  scale_linewidth_manual(name   = "Confidence Level", 
                         values = c(2, 1, 0.3)) +
  theme(legend.position   = "bottom",
        legend.margin     = margin(0, 50, 0, -20),
        legend.box.margin = margin(-5, -10, -5, -10),
        legend.text       = element_text(face = 3),
        legend.title      = element_text(face = 3)) +
  guides(colour = guide_legend(title.position = "top", title.hjust = 0.5),
         linewidth   = guide_legend(title.position = "top", title.hjust = 0.5, 
                                    override.aes = list(colour = "#1377E2")))

wet_mvmt_plot

Warning: Removed 5 rows containing missing values or values outside the scale range
(`geom_text()`).

ggsave(paste0("outputs/plots/pheno_model_wet_pooled_mvmt_", Sys.Date(), ".png"), device = "png", 
       width = 150, height = 150, units = "mm", scale = 1, dpi = 600)

Warning: Removed 5 rows containing missing values or values outside the scale range
(`geom_text()`).

Drop the filter argument when plotting all covariates

wet_mvmt_plot <- ggplot(data = coef_df %>% filter(str_detect(model_covariate, "pheno_end")), 
                        aes(y = Covariate, x = Estimate, col = factor(Preference))) +
  geom_vline(xintercept = 0, color = "black", lty = 2, lwd = 0.3) +
  geom_point(shape = 3, size = 2.5) +
  geom_errorbarh(data = coef_df2_no_mvmt,
                 aes(xmin = LCI, xmax = UCI, linewidth = factor(Level)), 
                 height = 0, alpha  = 0.5) +
  geom_text(aes(label = Significance, hjust = 0.5, vjust = 0), 
            show.legend = F) +
  scale_y_discrete(limits = rev(order_no_mvmt)) + 
  theme_classic() +
  scale_x_continuous(limits = c(min_x, max_x), breaks = seq(-2.5, max_x, 0.5)) +
  coord_capped_cart(left = "both", bottom = "both") +
  labs(x = expression(beta*"-Estimate")) +
  scale_color_manual(name   = "Preference", 
                     values = c("#238DFF", "#004691")) +
  scale_linewidth_manual(name   = "Confidence Level", 
                         values = c(2, 1, 0.3)) +
  theme(legend.position   = "bottom",
        legend.margin     = margin(0, 50, 0, -20),
        legend.box.margin = margin(-5, -10, -5, -10),
        legend.text       = element_text(face = 3),
        legend.title      = element_text(face = 3)) +
  guides(colour = guide_legend(title.position = "top", title.hjust = 0.5),
         linewidth   = guide_legend(title.position = "top", title.hjust = 0.5, 
                                    override.aes = list(colour = "#1377E2")))

wet_mvmt_plot

ggsave(paste0("outputs/plots/pheno_model_wet_pooled_mvmt_WITHOUT_MVMT_", Sys.Date(), ".png"), device = "png", 
       width = 150, height = 60, units = "mm", scale = 1, dpi = 600)

Updating movement parameters, with vegetation interactions

# from the step generation script
tentative_shape <- 0.438167
tentative_scale <- 534.3507
tentative_kappa <- 0.1848126

sl_extent <- 1000 # metres
sl_dry <- data.frame(x = seq(0, sl_extent, length.out = 1001))

Dry season

coef_df_dry_mvmt_mvparams <- data.frame("covariate" = names(fixef(buffalo.pheno.dry.mvmt)$cond),
                           "scaled_coefficient" = coef(summary(buffalo.pheno.dry.mvmt))$cond[, "Estimate"])

# step length coefficients
sl_dry_mvmt_mvparams <- coef_df_dry_mvmt_mvparams %>% filter(grepl("sl_scaled", covariate) & !grepl("log", covariate)) 
sl_dry_mvmt_mvparams <- sl_dry_mvmt_mvparams %>% mutate(
  scale_factor = sl_scale_sd,
  nat_coef = scaled_coefficient / scale_factor,
  pheno = c("Intercept", "Floodplain", "DryGrassland", "OpenForest", "TreeGrass"))
  
# log step length coefficients
log_sl_dry_mvmt_mvparams <- coef_df_dry_mvmt_mvparams %>% filter(grepl("log_sl_scaled", covariate)) 
log_sl_dry_mvmt_mvparams <- log_sl_dry_mvmt_mvparams %>% mutate(
  scale_factor = log_sl_scale_sd,
  nat_coef = scaled_coefficient / scale_factor,
  pheno = c("Intercept", "Floodplain", "DryGrassland", "OpenForest", "TreeGrass"))

# turning angle coefficients
cos_ta_dry_mvmt_mvparams <- coef_df_dry_mvmt_mvparams %>% filter(grepl("cos_ta_scaled", covariate)) 
cos_ta_dry_mvmt_mvparams <- cos_ta_dry_mvmt_mvparams %>% mutate(
  scale_factor = cos_ta_scale_sd,
  nat_coef = scaled_coefficient / scale_factor,
  pheno = c("Intercept", "Floodplain", "DryGrassland", "OpenForest", "TreeGrass"))

tentative_mean_sl <- tentative_shape * tentative_scale

# shrub_grass
shrub_grass_mvmt_dry <- data.frame(
  "pheno" = "ShrubGrass",
  "updated_shape" = tentative_shape + log_sl_dry_mvmt_mvparams$nat_coef[1],
  "updated_scale" = 1/((1/tentative_scale) - sl_dry_mvmt_mvparams$nat_coef[1]),
  "updated_kappa" = tentative_kappa + cos_ta_dry_mvmt_mvparams$nat_coef[1]
)

shrub_grass_mean_sl_dry <- shrub_grass_mvmt_dry$updated_shape * shrub_grass_mvmt_dry$updated_scale

# floodplain
floodplain_mvmt_dry <- data.frame(
  "pheno" = "Floodplain",
  "updated_shape" = tentative_shape + log_sl_dry_mvmt_mvparams$nat_coef[1] + log_sl_dry_mvmt_mvparams$nat_coef[2],
  "updated_scale" = 1/((1/tentative_scale) - (sl_dry_mvmt_mvparams$nat_coef[1] + sl_dry_mvmt_mvparams$nat_coef[2])),
  "updated_kappa" = tentative_kappa + cos_ta_dry_mvmt_mvparams$nat_coef[1] + cos_ta_dry_mvmt_mvparams$nat_coef[2]
)

floodplain_mean_sl_dry <- floodplain_mvmt_dry$updated_shape * floodplain_mvmt_dry$updated_scale

# dry_grassland
dry_grassland_mvmt_dry <- data.frame(
  "pheno" = "DryGrassland",
  "updated_shape" = tentative_shape + log_sl_dry_mvmt_mvparams$nat_coef[1] + log_sl_dry_mvmt_mvparams$nat_coef[3],
  "updated_scale" = 1/((1/tentative_scale) - (sl_dry_mvmt_mvparams$nat_coef[1] + sl_dry_mvmt_mvparams$nat_coef[3])),
  "updated_kappa" = tentative_kappa + cos_ta_dry_mvmt_mvparams$nat_coef[1] + cos_ta_dry_mvmt_mvparams$nat_coef[3]
)

dry_grassland_mean_sl_dry <- dry_grassland_mvmt_dry$updated_shape * dry_grassland_mvmt_dry$updated_scale

# open_forest
open_forest_mvmt_dry <- data.frame(
  "pheno" = "OpenForest",
  "updated_shape" = tentative_shape + log_sl_dry_mvmt_mvparams$nat_coef[1] + log_sl_dry_mvmt_mvparams$nat_coef[4],
  "updated_scale" = 1/((1/tentative_scale) - (sl_dry_mvmt_mvparams$nat_coef[1] + sl_dry_mvmt_mvparams$nat_coef[4])),
  "updated_kappa" = tentative_kappa + cos_ta_dry_mvmt_mvparams$nat_coef[1] + cos_ta_dry_mvmt_mvparams$nat_coef[4]
)

open_forest_mean_sl_dry <- open_forest_mvmt_dry$updated_shape * open_forest_mvmt_dry$updated_scale

# tree_grass
tree_grass_mvmt_dry <- data.frame(
  "pheno" = "TreeGrass",
  "updated_shape" = tentative_shape + log_sl_dry_mvmt_mvparams$nat_coef[1] + log_sl_dry_mvmt_mvparams$nat_coef[5],
  "updated_scale" = 1/((1/tentative_scale) - (sl_dry_mvmt_mvparams$nat_coef[1] + sl_dry_mvmt_mvparams$nat_coef[5])),
  "updated_kappa" = tentative_kappa + cos_ta_dry_mvmt_mvparams$nat_coef[1] + cos_ta_dry_mvmt_mvparams$nat_coef[5]
)

tree_grass_mean_sl_dry <- tree_grass_mvmt_dry$updated_shape * tree_grass_mvmt_dry$updated_scale

Creating Gamma distributions

sl_dry_mvmt <- data.frame(x = seq(0, sl_extent, length.out = 1001))

# Tentative distribution
sl_dry_mvmt$Tentative <- dgamma(x = sl_dry_mvmt$x,
                                 shape = tentative_shape,
                                 scale = tentative_scale)
# Shrub grass
sl_dry_mvmt$ShrubGrass <- dgamma(x = sl_dry_mvmt$x,
                             shape = shrub_grass_mvmt_dry$updated_shape,
                             scale = shrub_grass_mvmt_dry$updated_scale)
# Bare ground
sl_dry_mvmt$Floodplain <- dgamma(x = sl_dry_mvmt$x,
                             shape = floodplain_mvmt_dry$updated_shape,
                             scale = floodplain_mvmt_dry$updated_scale)
# Dry grassland
sl_dry_mvmt$DryGrassland <- dgamma(x = sl_dry_mvmt$x,
                             shape = dry_grassland_mvmt_dry$updated_shape,
                             scale = dry_grassland_mvmt_dry$updated_scale)
# Open forest
sl_dry_mvmt$OpenForest <- dgamma(x = sl_dry_mvmt$x,
                             shape = open_forest_mvmt_dry$updated_shape,
                             scale = open_forest_mvmt_dry$updated_scale)
# Tree grass
sl_dry_mvmt$TreeGrass <- dgamma(x = sl_dry_mvmt$x,
                             shape = tree_grass_mvmt_dry$updated_shape,
                             scale = tree_grass_mvmt_dry$updated_scale)

# Pivot from wide to long data
sl_dry_mvmt <- sl_dry_mvmt %>% 
  pivot_longer(cols = -x)

# Specify the order in which the coefficients should be plotted
order <- c(
  "ShrubGrass",
  "Floodplain",
  "DryGrassland",
  "OpenForest",
  "TreeGrass",
  "Tentative"
)

sl_dry_mvmt$name <- factor(sl_dry_mvmt$name, levels = order)


mean_values <- c("ShrubGrass" = round(shrub_grass_mean_sl_dry, 1),
                 "Floodplain" = round(floodplain_mean_sl_dry, 1), 
                 "DryGrassland" = round(dry_grassland_mean_sl_dry, 1), 
                 "OpenForest" = round(open_forest_mean_sl_dry, 1), 
                 "TreeGrass" = round(tree_grass_mean_sl_dry, 1), 
                 "Tentative" = round(tentative_mean_sl, 1))

# Construct dynamic labels
labels <- paste(names(mean_values), ": Mean SL = ", mean_values, sep = "")

# Plot
ggplot(sl_dry_mvmt, aes(x = x, y = value, color = name, linetype = name)) +
  geom_line(alpha = 0.75, linewidth = 1) +
  scale_x_continuous("Step Length (m)", limits = c(0,1000)) +
  scale_y_continuous("Probability Density", limits = c(0, 0.01)) +
  scale_colour_manual("Habitat",
                      values = c("ShrubGrass" = "#999900",
                                 "Floodplain" = "#3399FF",
                                 "DryGrassland" = "#CC9900",
                                 "OpenForest" = "#006600",
                                 "TreeGrass" = "#00CC66",
                                 "Tentative" = "black"),
                      labels = labels) +
  scale_linetype_manual("Habitat", 
                        values = c("ShrubGrass" = "solid", 
                                   "Floodplain" = "solid", 
                                   "DryGrassland" = "solid", 
                                   "OpenForest" = "solid", 
                                   "TreeGrass" = "solid", 
                                   "Tentative" = "dotted"),
                        labels = labels) +
  theme_bw() +
  theme(legend.position = c(0.65, 0.65))

Warning: A numeric `legend.position` argument in `theme()` was deprecated in ggplot2
3.5.0.
ℹ Please use the `legend.position.inside` argument of `theme()` instead.

ggsave(paste0("outputs/plots/pheno_model_sl_dry_pooled_mvmt_", Sys.Date(), ".png"), device = "png", 
       width = 150, height = 100, units = "mm", scale = 1, dpi = 1000)

Creating von Mises distributions

ta_dry_mvmt <- data.frame(x = seq(-pi, pi, length.out = 101))

# Tentative distribution
ta_dry_mvmt$Tentative <- circular::dvonmises(x = ta_dry_mvmt$x,
                                          mu = 0,
                                          kappa = tentative_kappa)
# Shrub grass
ta_dry_mvmt$ShrubGrass <- circular::dvonmises(x = ta_dry_mvmt$x,
                                          mu = 0,
                                          kappa = shrub_grass_mvmt_dry$updated_kappa)
# Bare ground
ta_dry_mvmt$Floodplain <- circular::dvonmises(x = ta_dry_mvmt$x,
                                          mu = 0,
                                          kappa = floodplain_mvmt_dry$updated_kappa)
# Dry grassland
ta_dry_mvmt$DryGrassland <- circular::dvonmises(x = ta_dry_mvmt$x,
                                          mu = 0,
                                          kappa = dry_grassland_mvmt_dry$updated_kappa)
# Open forest
ta_dry_mvmt$OpenForest <- circular::dvonmises(x = ta_dry_mvmt$x,
                                          mu = 0,
                                          kappa = open_forest_mvmt_dry$updated_kappa)
# Tree grass
ta_dry_mvmt$TreeGrass <- circular::dvonmises(x = ta_dry_mvmt$x,
                                          mu = 0,
                                          kappa = tree_grass_mvmt_dry$updated_kappa)

# Pivot from wide to long data
ta_dry_mvmt <- ta_dry_mvmt %>% 
  pivot_longer(cols = -x)

ta_dry_mvmt$name <- factor(ta_dry_mvmt$name, levels = order)

# Plot
ggplot(ta_dry_mvmt, 
       aes(x = x, y = value, color = factor(name), linetype = factor(name))
       ) +
  geom_line(alpha = 0.75, linewidth = 1) +
  xlab("Turning angle (radians)") +
  scale_y_continuous("Probability Density", limits = c(0.10, 0.215)) +
  scale_x_continuous(breaks = c(-pi, -pi/2, 0, pi/2, pi),
                     labels = c(expression(-pi, -pi/2, 0, pi/2, pi))) +
  scale_colour_manual("Habitat",
                      values = c("ShrubGrass" = "#999900",
                                 "Floodplain" = "#3399FF",
                                 "DryGrassland" = "#CC9900",
                                 "OpenForest" = "#006600",
                                 "TreeGrass" = "#00CC66",
                                 "Tentative" = "black")) +
  scale_linetype_manual("Habitat", 
                        values = c("ShrubGrass" = "solid", 
                                   "Floodplain" = "solid", 
                                   "DryGrassland" = "solid", 
                                   "OpenForest" = "solid", 
                                   "TreeGrass" = "solid", 
                                   "Tentative" = "dotted")) +
  # scale_colour_viridis_d("Habitat", option = "D") +
  # ggtitle("Dry season") +
  theme_bw() +
  theme(legend.position = "bottom")

ggsave(paste0("outputs/plots/pheno_model_ta_dry_pooled_mvmt_", Sys.Date(), ".png"), device = "png", 
       width = 150, height = 100, units = "mm", scale = 1, dpi = 600)

Wet season

coef_df_wet_mvmt_mvparams <- data.frame("covariate" = names(fixef(buffalo.pheno.wet.mvmt)$cond),
                                        "scaled_coefficient" = coef(summary(buffalo.pheno.wet.mvmt))$cond[, "Estimate"])

# step length coefficients
sl_wet_mvmt_mvparams <- coef_df_wet_mvmt_mvparams %>% filter(grepl("sl_scaled", covariate) & !grepl("log", covariate)) 
sl_wet_mvmt_mvparams <- sl_wet_mvmt_mvparams %>% mutate(
  scale_factor = sl_scale_sd,
  nat_coef = scaled_coefficient / scale_factor,
  pheno = c("Intercept", "Floodplain", "DryGrassland", "OpenForest", "TreeGrass"))

# log step length coefficients
log_sl_wet_mvmt_mvparams <- coef_df_wet_mvmt_mvparams %>% filter(grepl("log_sl_scaled", covariate)) 
log_sl_wet_mvmt_mvparams <- log_sl_wet_mvmt_mvparams %>% mutate(
  scale_factor = log_sl_scale_sd,
  nat_coef = scaled_coefficient / scale_factor,
  pheno = c("Intercept", "Floodplain", "DryGrassland", "OpenForest", "TreeGrass"))

# turning angle coefficients
cos_ta_wet_mvmt_mvparams <- coef_df_wet_mvmt_mvparams %>% filter(grepl("cos_ta_scaled", covariate)) 
cos_ta_wet_mvmt_mvparams <- cos_ta_wet_mvmt_mvparams %>% mutate(
  scale_factor = cos_ta_scale_sd,
  nat_coef = scaled_coefficient / scale_factor,
  pheno = c("Intercept", "Floodplain", "DryGrassland", "OpenForest", "TreeGrass"))

# shrub_grass
shrub_grass_mvmt_wet <- data.frame(
  "pheno" = "ShrubGrass",
  "updated_shape" = tentative_shape + log_sl_wet_mvmt_mvparams$nat_coef[1],
  "updated_scale" = 1/((1/tentative_scale) - sl_wet_mvmt_mvparams$nat_coef[1]),
  "updated_kappa" = tentative_kappa + cos_ta_wet_mvmt_mvparams$nat_coef[1]
)

shrub_grass_mean_sl_wet <- shrub_grass_mvmt_wet$updated_shape * shrub_grass_mvmt_wet$updated_scale

# floodplain
floodplain_mvmt_wet <- data.frame(
  "pheno" = "Floodplain",
  "updated_shape" = tentative_shape + log_sl_wet_mvmt_mvparams$nat_coef[1] + log_sl_wet_mvmt_mvparams$nat_coef[2],
  "updated_scale" = 1/((1/tentative_scale) - (sl_wet_mvmt_mvparams$nat_coef[1] + sl_wet_mvmt_mvparams$nat_coef[2])),
  "updated_kappa" = tentative_kappa + cos_ta_wet_mvmt_mvparams$nat_coef[1] + cos_ta_wet_mvmt_mvparams$nat_coef[2]
)

floodplain_mean_sl_wet <- floodplain_mvmt_wet$updated_shape * floodplain_mvmt_wet$updated_scale

# dry_grassland
dry_grassland_mvmt_wet <- data.frame(
  "pheno" = "wetGrassland",
  "updated_shape" = tentative_shape + log_sl_wet_mvmt_mvparams$nat_coef[1] + log_sl_wet_mvmt_mvparams$nat_coef[3],
  "updated_scale" = 1/((1/tentative_scale) - (sl_wet_mvmt_mvparams$nat_coef[1] + sl_wet_mvmt_mvparams$nat_coef[3])),
  "updated_kappa" = tentative_kappa + cos_ta_wet_mvmt_mvparams$nat_coef[1] + cos_ta_wet_mvmt_mvparams$nat_coef[3]
)

dry_grassland_mean_sl_wet <- dry_grassland_mvmt_wet$updated_shape * dry_grassland_mvmt_wet$updated_scale

# open_forest
open_forest_mvmt_wet <- data.frame(
  "pheno" = "OpenForest",
  "updated_shape" = tentative_shape + log_sl_wet_mvmt_mvparams$nat_coef[1] + log_sl_wet_mvmt_mvparams$nat_coef[4],
  "updated_scale" = 1/((1/tentative_scale) - (sl_wet_mvmt_mvparams$nat_coef[1] + sl_wet_mvmt_mvparams$nat_coef[4])),
  "updated_kappa" = tentative_kappa + cos_ta_wet_mvmt_mvparams$nat_coef[1] + cos_ta_wet_mvmt_mvparams$nat_coef[4]
)

open_forest_mean_sl_wet <- open_forest_mvmt_wet$updated_shape * open_forest_mvmt_wet$updated_scale

# tree_grass
tree_grass_mvmt_wet <- data.frame(
  "pheno" = "TreeGrass",
  "updated_shape" = tentative_shape + log_sl_wet_mvmt_mvparams$nat_coef[1] + log_sl_wet_mvmt_mvparams$nat_coef[5],
  "updated_scale" = 1/((1/tentative_scale) - (sl_wet_mvmt_mvparams$nat_coef[1] + sl_wet_mvmt_mvparams$nat_coef[5])),
  "updated_kappa" = tentative_kappa + cos_ta_wet_mvmt_mvparams$nat_coef[1] + cos_ta_wet_mvmt_mvparams$nat_coef[5]
)

tree_grass_mean_sl_wet <- tree_grass_mvmt_wet$updated_shape * tree_grass_mvmt_wet$updated_scale

Creating Gamma distributions

sl_wet_mvmt <- data.frame(x = seq(0, sl_extent, length.out = 1001))

# Tentative distribution
sl_wet_mvmt$Tentative <- dgamma(x = sl_wet_mvmt$x,
                                shape = tentative_shape,
                                scale = tentative_scale)
# Shrub grass
sl_wet_mvmt$ShrubGrass <- dgamma(x = sl_wet_mvmt$x,
                                 shape = shrub_grass_mvmt_wet$updated_shape,
                                 scale = shrub_grass_mvmt_wet$updated_scale)
# Bare ground
sl_wet_mvmt$Floodplain <- dgamma(x = sl_wet_mvmt$x,
                                 shape = floodplain_mvmt_wet$updated_shape,
                                 scale = floodplain_mvmt_wet$updated_scale)
# Dry grassland
sl_wet_mvmt$DryGrassland <- dgamma(x = sl_wet_mvmt$x,
                                   shape = dry_grassland_mvmt_wet$updated_shape,
                                   scale = dry_grassland_mvmt_wet$updated_scale)
# Open forest
sl_wet_mvmt$OpenForest <- dgamma(x = sl_wet_mvmt$x,
                                 shape = open_forest_mvmt_wet$updated_shape,
                                 scale = open_forest_mvmt_wet$updated_scale)
# Tree grass
sl_wet_mvmt$TreeGrass <- dgamma(x = sl_wet_mvmt$x,
                                shape = tree_grass_mvmt_wet$updated_shape,
                                scale = tree_grass_mvmt_wet$updated_scale)

# Pivot from wide to long data
sl_wet_mvmt <- sl_wet_mvmt %>% 
  pivot_longer(cols = -x)

sl_wet_mvmt$name <- factor(sl_wet_mvmt$name, levels = order)

mean_values_wet <- c("ShrubGrass" = round(shrub_grass_mean_sl_wet, 1),
                     "Floodplain" = round(floodplain_mean_sl_wet, 1), 
                     "DryGrassland" = round(dry_grassland_mean_sl_wet, 1), 
                     "OpenForest" = round(open_forest_mean_sl_wet, 1), 
                     "TreeGrass" = round(tree_grass_mean_sl_wet, 1), 
                     "Tentative" = round(tentative_mean_sl, 1))

# Construct dynamic labels
labels <- paste(names(mean_values_wet), ": Mean SL = ", mean_values_wet, sep = "")

# Plot
ggplot(sl_wet_mvmt, aes(x = x, y = value, color = factor(name), linetype = name)) +
  geom_line(alpha = 0.75, linewidth = 1) +
  scale_x_continuous("Step Length (m)", limits = c(0,1000)) +
  scale_y_continuous("Probability Density", limits = c(0, 0.01)) +
  scale_colour_manual("Habitat",
                      values = c("ShrubGrass" = "#999900",
                                 "Floodplain" = "#3399FF",
                                 "DryGrassland" = "#CC9900",
                                 "OpenForest" = "#006600",
                                 "TreeGrass" = "#00CC66",
                                 "Tentative" = "black"),
                      labels = labels) +
  scale_linetype_manual("Habitat", 
                        values = c("ShrubGrass" = "solid", 
                                   "Floodplain" = "solid", 
                                   "DryGrassland" = "solid", 
                                   "OpenForest" = "solid", 
                                   "TreeGrass" = "solid", 
                                   "Tentative" = "dotted"),
                        labels = labels) +
  # scale_colour_viridis_d("Habitat", option = "D") +
  # ggtitle("Wet season") +
  theme_bw() +
  theme(legend.position = c(0.65, 0.65))

ggsave(paste0("outputs/plots/pheno_model_sl_wet_pooled_mvmt_", Sys.Date(), ".png"), device = "png", 
       width = 150, height = 100, units = "mm", scale = 1, dpi = 1000)

Creating von Mises distributions

ta_wet_mvmt <- data.frame(x = seq(-pi, pi, length.out = 101))

# Tentative distribution
ta_wet_mvmt$Tentative <- circular::dvonmises(x = ta_wet_mvmt$x,
                                             mu = 0,
                                             kappa = tentative_kappa)
# Shrub grass
ta_wet_mvmt$ShrubGrass <- circular::dvonmises(x = ta_wet_mvmt$x,
                                              mu = 0,
                                              kappa = shrub_grass_mvmt_wet$updated_kappa)
# Bare ground
ta_wet_mvmt$Floodplain <- circular::dvonmises(x = ta_wet_mvmt$x,
                                              mu = 0,
                                              kappa = floodplain_mvmt_wet$updated_kappa)
# Dry grassland
ta_wet_mvmt$DryGrassland <- circular::dvonmises(x = ta_wet_mvmt$x,
                                                mu = 0,
                                                kappa = dry_grassland_mvmt_wet$updated_kappa)
# Open forest
ta_wet_mvmt$OpenForest <- circular::dvonmises(x = ta_wet_mvmt$x,
                                              mu = 0,
                                              kappa = open_forest_mvmt_wet$updated_kappa)
# Tree grass
ta_wet_mvmt$TreeGrass <- circular::dvonmises(x = ta_wet_mvmt$x,
                                             mu = 0,
                                             kappa = tree_grass_mvmt_wet$updated_kappa)

# Pivot from wide to long data
ta_wet_mvmt <- ta_wet_mvmt %>% 
  pivot_longer(cols = -x)

ta_wet_mvmt$name <- factor(ta_wet_mvmt$name, levels = order)

# Plot
ggplot(ta_wet_mvmt, aes(x = x, y = value, color = factor(name), linetype = name)) +
  geom_line(alpha = 0.75, linewidth = 1) +
  xlab("Turning angle (radians)") +
  scale_y_continuous("Probability Density", limits = c(0.10, 0.215)) +
  scale_x_continuous(breaks = c(-pi, -pi/2, 0, pi/2, pi),
                     labels = c(expression(-pi, -pi/2, 0, pi/2, pi))) +
  scale_colour_manual("Habitat",
                      values = c("ShrubGrass" = "#999900",
                                 "Floodplain" = "#3399FF",
                                 "DryGrassland" = "#CC9900",
                                 "OpenForest" = "#006600",
                                 "TreeGrass" = "#00CC66",
                                 "Tentative" = "black")) +
  scale_linetype_manual("Habitat", 
                        values = c("ShrubGrass" = "solid", 
                                   "Floodplain" = "solid", 
                                   "DryGrassland" = "solid", 
                                   "OpenForest" = "solid", 
                                   "TreeGrass" = "solid", 
                                   "Tentative" = "dotted")) +
  # scale_colour_viridis_d("Habitat", option = "D") +
  # ggtitle("Dry season") +
  theme_bw() +
  theme(legend.position = "bottom")

ggsave(paste0("outputs/plots/pheno_model_ta_wet_pooled_mvmt_", Sys.Date(), ".png"), device = "png", 
       width = 150, height = 100, units = "mm", scale = 1, dpi = 600)

Session Information

sessionInfo()

R version 4.5.0 (2025-04-11)
Platform: aarch64-apple-darwin20
Running under: macOS Sequoia 15.4.1

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRblas.0.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

time zone: Australia/Brisbane
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] amt_0.2.2.0         circular_0.5-1      lemon_0.5.0        
 [4] beepr_2.0           tictoc_1.2.1        glmmTMB_1.1.11     
 [7] sjPlot_2.8.17       TwoStepCLogit_1.2.6 lubridate_1.9.4    
[10] forcats_1.0.0       stringr_1.5.1       dplyr_1.1.4        
[13] purrr_1.0.4         readr_2.1.5         tidyr_1.3.1        
[16] tibble_3.2.1        ggplot2_3.5.2       tidyverse_2.0.0    

loaded via a namespace (and not attached):
 [1] sjlabelled_1.2.0    tidyselect_1.2.1    farver_2.1.2       
 [4] fastmap_1.2.0       bayestestR_0.15.3   sjstats_0.19.0     
 [7] digest_0.6.37       timechange_0.3.0    lifecycle_1.0.4    
[10] sf_1.0-20           survival_3.8-3      magrittr_2.0.3     
[13] compiler_4.5.0      rlang_1.1.6         tools_4.5.0        
[16] yaml_2.3.10         knitr_1.50          labeling_0.4.3     
[19] bit_4.6.0           classInt_0.4-11     plyr_1.8.9         
[22] RColorBrewer_1.1-3  KernSmooth_2.23-26  withr_3.0.2        
[25] numDeriv_2016.8-1.1 grid_4.5.0          datawizard_1.1.0   
[28] e1071_1.7-16        scales_1.4.0        MASS_7.3-65        
[31] insight_1.2.0       cli_3.6.5           mvtnorm_1.3-3      
[34] crayon_1.5.3        rmarkdown_2.29      ragg_1.4.0         
[37] reformulas_0.4.1    generics_0.1.3      rstudioapi_0.17.1  
[40] performance_0.13.0  tzdb_0.5.0          parameters_0.25.0  
[43] minqa_1.2.8         DBI_1.2.3           proxy_0.4-27       
[46] audio_0.1-11        splines_4.5.0       parallel_4.5.0     
[49] effectsize_1.0.0    vctrs_0.6.5         boot_1.3-31        
[52] Matrix_1.7-3        jsonlite_2.0.0      hms_1.1.3          
[55] bit64_4.6.0-1       systemfonts_1.2.3   units_0.8-7        
[58] glue_1.8.0          nloptr_2.2.1        stringi_1.8.7      
[61] gtable_0.3.6        ggeffects_2.2.1     lme4_1.1-37        
[64] pillar_1.10.2       htmltools_0.5.8.1   R6_2.6.1           
[67] TMB_1.9.17          textshaping_1.0.1   Rdpack_2.6.4       
[70] vroom_1.6.5         evaluate_1.0.3      lattice_0.22-6     
[73] rbibutils_2.3       class_7.3-23        Rcpp_1.0.14        
[76] gridExtra_2.3       nlme_3.1-168        mgcv_1.9-1         
[79] xfun_0.52           sjmisc_2.8.10       pkgconfig_2.0.3